library(readxl)
library(dplyr)
library(glmtoolbox) # Para stepCriterion() y envelope()
library(aplore3) # Para el dataset burn1000
library(ISLR2) # Para el dataset Smarket15 Modelos Lineales Generalizados GAM
16 Introducción
Este documento explora los Modelos Lineales Generalizados (GLM), con un enfoque en la regresión logística para variables de respuesta binarias. Se divide en tres ejercicios principales:
Ejercicio 1: Se ajusta un modelo logístico, se interpretan sus coeficientes, se comparan diferentes funciones de enlace y se realiza una selección automática de variables.
Ejercicio 2: Se aplica la regresión logística para predecir la dirección del mercado de valores, introduciendo el concepto de dividir los datos en conjuntos de entrenamiento y prueba.
Ejercicio 3: Se profundiza en la evaluación de la habilidad predictiva del modelo, utilizando muestreo estratificado, curvas ROC, y la optimización de un punto de corte (
tau) para la clasificación.
16.0.1 Librerías Requeridas
17 EJERCICIO 1: Bondad de Ajuste en Regresión Logística
Se utiliza el conjunto de datos burn1000 para modelar la probabilidad de muerte (death) en función de varios factores de riesgo.
17.1 a) Análisis de Estadística Descriptiva
burn1000 <- aplore3::burn1000
# Se crea la variable de respuesta binaria: 1 = "Dead" (éxito), 0 = "Alive"
burn1000 <- within(burn1000, death2 <- ifelse(death == "Dead", 1, 0))
str(burn1000)'data.frame': 1000 obs. of 10 variables:
$ id : int 1 2 3 4 5 6 7 8 9 10 ...
$ facility: int 11 1 12 1 1 6 22 1 1 1 ...
$ death : Factor w/ 2 levels "Alive","Dead": 1 1 1 1 1 1 1 1 1 1 ...
$ age : num 26.6 2 22 37.3 52.1 50.2 2.5 53.8 31.9 41.1 ...
$ gender : Factor w/ 2 levels "Female","Male": 2 1 1 2 2 2 1 1 2 2 ...
$ race : Factor w/ 2 levels "Non-White","White": 2 1 1 2 2 2 1 2 2 2 ...
$ tbsa : num 25.3 5 2 2 6 7 7 0.9 2 22 ...
$ inh_inj : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
$ flame : Factor w/ 2 levels "No","Yes": 2 1 1 1 2 1 1 2 1 2 ...
$ death2 : num 0 0 0 0 0 0 0 0 0 0 ...summary(burn1000) id facility death age gender
Min. : 1.0 Min. : 1.00 Alive:850 Min. : 0.10 Female:295
1st Qu.: 250.8 1st Qu.: 2.00 Dead :150 1st Qu.:10.85 Male :705
Median : 500.5 Median : 8.00 Median :31.95
Mean : 500.5 Mean :11.56 Mean :33.29
3rd Qu.: 750.2 3rd Qu.:18.25 3rd Qu.:51.23
Max. :1000.0 Max. :40.00 Max. :89.70
race tbsa inh_inj flame death2
Non-White:411 Min. : 0.10 No :878 No :471 Min. :0.00
White :589 1st Qu.: 2.50 Yes:122 Yes:529 1st Qu.:0.00
Median : 6.00 Median :0.00
Mean :13.54 Mean :0.15
3rd Qu.:16.00 3rd Qu.:0.00
Max. :98.00 Max. :1.00 17.2 b-d) Ajuste de un Primer Modelo Logístico
Se ajusta un GLM con una distribución binomial y una función de enlace logit. El modelo es: \[
\log\left(\frac{p_i}{1-p_i}\right) = \beta_0 + \beta_1 x_{i1} + \dots + \beta_k x_{ik}
\] donde \(p_i\) es la probabilidad de “éxito” (muerte) para el individuo \(i\).
# Definir el Modelo Lineal Generalizado
fit1 <- glm(death2 ~ age + gender + race + tbsa + inh_inj + flame,
family = binomial("logit"), data = burn1000)
summary(fit1)
Call:
glm(formula = death2 ~ age + gender + race + tbsa + inh_inj +
flame, family = binomial("logit"), data = burn1000)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -7.695153 0.691169 -11.134 < 2e-16 ***
age 0.082890 0.008629 9.606 < 2e-16 ***
genderMale -0.201494 0.307784 -0.655 0.512687
raceWhite -0.701389 0.309781 -2.264 0.023565 *
tbsa 0.089345 0.009087 9.832 < 2e-16 ***
inh_injYes 1.365277 0.361780 3.774 0.000161 ***
flameYes 0.582578 0.354493 1.643 0.100298
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 845.42 on 999 degrees of freedom
Residual deviance: 336.46 on 993 degrees of freedom
AIC: 350.46
Number of Fisher Scoring iterations: 717.3 e) Interpretación de los Parámetros (Odds Ratios)
En un modelo logístico, los coeficientes exponenciados, \(e^{\beta_j}\), se interpretan como Odds Ratios (razón de chances).
# Odds Ratios puntuales
exp(coef(fit1)[-1]) age genderMale raceWhite tbsa inh_injYes flameYes
1.0864226 0.8175088 0.4958961 1.0934576 3.9168097 1.7906491 # Un OR > 1 indica que el chance de éxito aumenta.
# Un OR < 1 indica que el chance de éxito disminuye.
# Cambio porcentual en los odds: (OR - 1) * 100
(exp(coef(fit1)[-1]) - 1) * 100 age genderMale raceWhite tbsa inh_injYes flameYes
8.642262 -18.249119 -50.410395 9.345756 291.680974 79.064907 # Intervalos de confianza del 95% para los Odds Ratios
exp(confint(fit1)[-1, ])Waiting for profiling to be done... 2.5 % 97.5 %
age 1.0691824 1.1060770
genderMale 0.4486691 1.5047454
raceWhite 0.2679382 0.9063323
tbsa 1.0752625 1.1143712
inh_injYes 1.9331829 8.0230447
flameYes 0.9052347 3.657049017.4 f) Elección de una Función de Enlace por Bondad de Ajuste
Se comparan cuatro funciones de enlace para la distribución binomial (logit, probit, cloglog, cauchit) utilizando los criterios de información AIC y BIC. Valores más bajos indican un mejor ajuste.
fit2 <- update(fit1, family = binomial("probit"))
fit3 <- update(fit1, family = binomial("cloglog"))Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredfit4 <- update(fit1, family = binomial("cauchit"))
# Comparación de modelos
cbind(AIC(fit1, fit2, fit3, fit4), BIC(fit1, fit2, fit3, fit4)) df AIC df BIC
fit1 7 350.4623 7 384.8166
fit2 7 347.0558 7 381.4101
fit3 7 374.6395 7 408.9938
fit4 7 389.5935 7 423.9478Conclusión: Según AIC y BIC, el modelo fit2 (con enlace probit) ofrece la mejor bondad de ajuste.
17.5 g-h) Selección Automática del Modelo
Se utiliza el método de eliminación hacia atrás (backward) basado en el criterio BIC para encontrar el subconjunto de predictores más parsimonioso, partiendo de un modelo con interacciones.
fit5 <- glm(death2 ~ age + gender + race + tbsa + inh_inj + flame + age*inh_inj + tbsa*inh_inj,
family = binomial("logit"), data = burn1000)
a <- stepCriterion(fit5, direction = "backward", criterion = "bic")
Family: binomial
Link function: logit
Criterion: BIC
Initial model:
~ age + gender + race + tbsa + inh_inj + flame + age * inh_inj + tbsa * inh_inj
Step 0 :
df AIC BIC adj.R-squared P(Chisq>)(*)
- gender 1 327.92 367.18 0.6284 0.933130
- flame 1 330.74 370.00 0.6251 0.097426
- race 1 334.78 374.04 0.6203 0.009797
<none> 329.91 374.08 0.6281
- inh_inj:tbsa 1 337.59 376.85 0.6169 0.002061
- age:inh_inj 1 351.79 391.06 0.6000 2.31e-06
Step 1 : - gender
df AIC BIC adj.R-squared P(Chisq>)(*)
- flame 1 328.74 363.09 0.6255 0.097813
<none> 327.92 367.18 0.6284
- race 1 332.83 367.19 0.6206 0.009622
- inh_inj:tbsa 1 335.61 369.96 0.6173 0.002051
- age:inh_inj 1 350.36 384.71 0.5997 1.643e-06
Step 2 : - flame
df AIC BIC adj.R-squared P(Chisq>)(*)
- race 1 332.42 361.87 0.6191 0.01843
<none> 328.74 363.09 0.6255
- inh_inj:tbsa 1 337.19 366.64 0.6134 0.00142
+ gender 1 330.74 370.00 0.6251 0.96626
- age:inh_inj 1 351.00 380.44 0.5970 1.744e-06
Step 3 : - race
df AIC BIC adj.R-squared P(Chisq>)(*)
<none> 332.42 361.87 0.6191
- inh_inj:tbsa 1 340.10 364.63 0.6080 0.002075
+ flame 1 332.83 367.19 0.6206 0.212331
+ gender 1 334.33 368.68 0.6188 0.758727
- age:inh_inj 1 353.46 378.00 0.5921 2.823e-06
Final model:
~ age + tbsa + inh_inj + age:inh_inj + tbsa:inh_inj
****************************************************************************
(*) p-values of the Wald testfit6 <- update(fit5, formula. = a$final)
summary(fit6)
Call:
glm(formula = death2 ~ age + tbsa + inh_inj + age:inh_inj + tbsa:inh_inj,
family = binomial("logit"), data = burn1000)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -10.31738 1.03278 -9.990 < 2e-16 ***
age 0.11561 0.01352 8.553 < 2e-16 ***
tbsa 0.11221 0.01439 7.797 6.36e-15 ***
inh_injYes 7.10938 1.25863 5.649 1.62e-08 ***
age:inh_injYes -0.08004 0.01709 -4.683 2.82e-06 ***
tbsa:inh_injYes -0.05538 0.01799 -3.079 0.00207 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 845.42 on 999 degrees of freedom
Residual deviance: 320.42 on 994 degrees of freedom
AIC: 332.42
Number of Fisher Scoring iterations: 8# Se re-evalúan las funciones de enlace con el nuevo modelo
fit7 <- update(fit6, family = binomial("probit"))
fit8 <- update(fit6, family = binomial("cloglog"))Warning: glm.fit: algorithm did not convergeWarning: glm.fit: fitted probabilities numerically 0 or 1 occurredfit9 <- update(fit6, family = binomial("cauchit"))
cbind(AIC(fit6, fit7, fit8, fit9), BIC(fit6, fit7, fit8, fit9)) df AIC df BIC
fit6 6 332.4233 6 361.8698
fit7 6 330.5678 6 360.0143
fit8 6 342.1889 6 371.6355
fit9 6 368.4801 6 397.9266Conclusión Final: El modelo fit7 (con predictores seleccionados y enlace probit) es el mejor según el criterio BIC.
17.6 i) Validación de Supuestos
Se utilizan gráficos de diagnóstico para evaluar el ajuste del modelo final.
# Gráfico de envelope para los residuales (similar a un Q-Q plot con bandas de confianza)
envelope(fit7)
|
| | 0%
|
|+ | 4%
|
|++ | 8%
|
|+++ | 12%
|
|++++ | 16%
|
|+++++ | 20%
|
|++++++ | 24%
|
|+++++++ | 28%
|
|++++++++ | 32%
|
|+++++++++ | 36%
|
|++++++++++ | 40%
|
|+++++++++++ | 44%
|
|++++++++++++ | 48%
|
|+++++++++++++ | 52%
|
|++++++++++++++ | 56%
|
|+++++++++++++++ | 60%
|
|++++++++++++++++ | 64%
|
|+++++++++++++++++ | 68%
|
|++++++++++++++++++ | 72%
|
|+++++++++++++++++++ | 76%
|
|++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++ | 84%
|
|++++++++++++++++++++++ | 88%
|
|+++++++++++++++++++++++ | 92%
|
|++++++++++++++++++++++++ | 96%
|
|+++++++++++++++++++++++++| 100%
18 EJERCICIO 2: Predicción del Mercado de Valores
Se utiliza el conjunto de datos Smarket para predecir la dirección del mercado (Up o Down).
18.1 b) Ajuste del Modelo
glm.fits <- glm(Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume,
data = Smarket, family = binomial)
summary(glm.fits)
Call:
glm(formula = Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 +
Volume, family = binomial, data = Smarket)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.126000 0.240736 -0.523 0.601
Lag1 -0.073074 0.050167 -1.457 0.145
Lag2 -0.042301 0.050086 -0.845 0.398
Lag3 0.011085 0.049939 0.222 0.824
Lag4 0.009359 0.049974 0.187 0.851
Lag5 0.010313 0.049511 0.208 0.835
Volume 0.135441 0.158360 0.855 0.392
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1731.2 on 1249 degrees of freedom
Residual deviance: 1727.6 on 1243 degrees of freedom
AIC: 1741.6
Number of Fisher Scoring iterations: 318.2 c) Matriz de Confusión (con todos los datos)
Se evalúa el rendimiento del modelo en los mismos datos con los que fue entrenado (esto generalmente sobreestima el rendimiento real).
glm.probs <- predict(glm.fits, type = "response")
glm.pred <- rep("Down", nrow(Smarket))
glm.pred[glm.probs > 0.5] <- "Up"
# Matriz de Confusión
table(Predicción = glm.pred, Real = Smarket$Direction) Real
Predicción Down Up
Down 145 141
Up 457 507# Tasa de acierto global (Accuracy)
mean(glm.pred == Smarket$Direction)[1] 0.521618.3 d) Entrenamiento y Prueba
Un enfoque más realista es dividir los datos: entrenar el modelo con datos antiguos (hasta 2004) y probar su rendimiento en datos nuevos (2005).
train <- (Smarket$Year < 2005)
Smarket.2005 <- Smarket[!train, ]
Direction.2005 <- Smarket$Direction[!train]
glm.fits.train <- glm(Direction ~ Lag1 + Lag2, # Modelo simplificado
data = Smarket, family = binomial, subset = train)
glm.probs.test <- predict(glm.fits.train, Smarket.2005, type = "response")
glm.pred.test <- rep("Down", nrow(Smarket.2005))
glm.pred.test[glm.probs.test > 0.5] <- "Up"
table(Predicción = glm.pred.test, Real = Direction.2005) Real
Predicción Down Up
Down 35 35
Up 76 106mean(glm.pred.test == Direction.2005)[1] 0.559523819 EJERCICIO 3: Habilidad Predictiva
Se vuelve al conjunto de datos burn1000 para realizar una evaluación predictiva más rigurosa.
19.1 a) Muestreo Estratificado
Dado que las clases (Dead/Alive) están desbalanceadas, se utiliza muestreo estratificado para asegurar que la proporción de clases sea la misma en los conjuntos de entrenamiento y prueba.
set.seed(123)
# Se toma un 70% de cada estrato para el conjunto de entrenamiento
stratified <- burn1000 %>%
group_by(death) %>%
slice_sample(prop = 0.7)
# El 30% restante forma el conjunto de prueba
test <- burn1000[setdiff(burn1000$id, stratified$id), ]
summary(burn1000$death)Alive Dead
850 150 summary(stratified$death)Alive Dead
595 105 19.2 b) Entrenamiento del Modelo y Selección de Variables
El proceso de ajuste y selección se realiza únicamente con el conjunto de entrenamiento.
# Se usa el modelo con interacciones como punto de partida
fit5b <- glm(death2 ~ age + gender + race + tbsa + inh_inj + flame + age*inh_inj + tbsa*inh_inj,
family = binomial("logit"), data = stratified)
a <- stepCriterion(fit5b, direction = "backward", criterion = "bic")
Family: binomial
Link function: logit
Criterion: BIC
Initial model:
~ age + gender + race + tbsa + inh_inj + flame + age * inh_inj + tbsa * inh_inj
Step 0 :
df AIC BIC adj.R-squared P(Chisq>)(*)
- gender 1 232.66 269.07 0.6302 0.892068
- flame 1 238.04 274.45 0.6210 0.023533
<none> 234.64 275.60 0.6297
- race 1 239.42 275.83 0.6186 0.010751
- inh_inj:tbsa 1 240.76 277.17 0.6164 0.005185
- age:inh_inj 1 253.72 290.12 0.5942 8.119e-06
Step 1 : - gender
df AIC BIC adj.R-squared P(Chisq>)(*)
- flame 1 236.04 267.90 0.6215 0.023837
<none> 232.66 269.07 0.6302
- race 1 237.46 269.32 0.6191 0.010725
- inh_inj:tbsa 1 238.79 270.64 0.6169 0.005162
- age:inh_inj 1 252.13 283.99 0.5941 6.321e-06
Step 2 : - flame
df AIC BIC adj.R-squared P(Chisq>)(*)
- race 1 238.68 265.98 0.6142 0.033530
<none> 236.04 267.90 0.6215
- inh_inj:tbsa 1 243.42 270.73 0.6061 0.002966
+ gender 1 238.04 274.45 0.6210 0.986257
- age:inh_inj 1 254.27 281.57 0.5877 9.978e-06
Step 3 : - race
df AIC BIC adj.R-squared P(Chisq>)(*)
<none> 238.68 265.98 0.6142
- inh_inj:tbsa 1 245.66 268.42 0.5995 0.00353
+ flame 1 237.46 269.32 0.6191 0.07887
+ gender 1 240.59 272.45 0.6138 0.76611
- age:inh_inj 1 255.86 278.62 0.5822 1.531e-05
Final model:
~ age + tbsa + inh_inj + age:inh_inj + tbsa:inh_inj
****************************************************************************
(*) p-values of the Wald testfit6b <- update(fit5b, formula. = a$final)
# Se comparan de nuevo las funciones de enlace
fit7b <- update(fit6b, family = binomial("probit"))
fit8b <- update(fit6b, family = binomial("cloglog"))Warning: glm.fit: algorithm did not convergeWarning: glm.fit: fitted probabilities numerically 0 or 1 occurredfit9b <- update(fit6b, family = binomial("cauchit"))19.3 c) Selección del Mejor Modelo Global por Habilidad Predictiva
Se utilizan las Curvas ROC (Receiver Operating Characteristic) para comparar el rendimiento predictivo de los modelos en el conjunto de prueba. El área bajo la curva (AUC) es una medida global del poder discriminatorio del modelo.
# Se obtienen las predicciones de probabilidad para el conjunto de prueba
pr6b <- predict(fit6b, newdata = test, type = "response")
pr7b <- predict(fit7b, newdata = test, type = "response")
pr8b <- predict(fit8b, newdata = test, type = "response")
pr9b <- predict(fit9b, newdata = test, type = "response")
# Se comparan las curvas ROC
testres <- test$death2
ROCc(cbind(testres, pr6b), main="Modelo 6 (logit)")
Area Under ROC Curve = 0.968, 95% CI: (0.949, 0.987)
Gini Coefficient = 0.936, 95% CI: (0.898, 0.974)
K-S Statistic = 0.834 ROCc(cbind(testres, pr7b), main="Modelo 7 (probit)")
Area Under ROC Curve = 0.968, 95% CI: (0.949, 0.987)
Gini Coefficient = 0.936, 95% CI: (0.898, 0.974)
K-S Statistic = 0.83 ROCc(cbind(testres, pr8b), main="Modelo 8 (cloglog)")
Area Under ROC Curve = 0.969, 95% CI: (0.95, 0.987)
Gini Coefficient = 0.937, 95% CI: (0.9, 0.974)
K-S Statistic = 0.816 ROCc(cbind(testres, pr9b), main="Modelo 9 (cauchit)")
Area Under ROC Curve = 0.966, 95% CI: (0.946, 0.986)
Gini Coefficient = 0.933, 95% CI: (0.893, 0.972)
K-S Statistic = 0.838 Conclusión: Se elige el modelo con el AUC más alto (en este caso, fit7b).
19.4 Elección del Punto de Corte (tau) Óptimo
Una vez elegido el mejor modelo (fit7b), se debe seleccionar un punto de corte (tau) para convertir las probabilidades predichas en clasificaciones (“Dead” / “Alive”). La elección de tau depende del objetivo (minimizar errores, maximizar precisión, etc.).
tau <- seq(0.1, 0.9, by = 0.05)
AER <- recall <- precision <- F1 <- NULL
exito <- "1"; frac <- "0"
for (i in 1:length(tau)){
glm.pred <- rep("0", length(testres))
glm.pred[pr7b > tau[i]] <- "1"
tab<-table(glm.pred, testres)
if (!frac %in% rownames(tab)){
tab<-rbind(tab,c(0,0))
rownames(tab)[2]<-frac
}
if (!exito %in% rownames(tab)){
tab<-rbind(tab,c(0,0))
rownames(tab)[2]<-exito
}
AER[i]<-(tab[exito,frac]+tab[frac,exito])/sum(tab)
precision[i]<-(tab[exito,exito])/sum(tab[exito,])
recall[i]<-(tab[exito,exito])/sum(tab[,exito])
F1[i]<-2*precision[i]*recall[i]/(precision[i]+recall[i])
}
cbind(tau,AER,precision,recall,F1) tau AER precision recall F1
[1,] 0.10 0.14000000 0.5180723 0.9555556 0.6718750
[2,] 0.15 0.11666667 0.5694444 0.9111111 0.7008547
[3,] 0.20 0.10333333 0.6060606 0.8888889 0.7207207
[4,] 0.25 0.07666667 0.7115385 0.8222222 0.7628866
[5,] 0.30 0.07666667 0.7200000 0.8000000 0.7578947
[6,] 0.35 0.07000000 0.7608696 0.7777778 0.7692308
[7,] 0.40 0.06333333 0.8250000 0.7333333 0.7764706
[8,] 0.45 0.06333333 0.8421053 0.7111111 0.7710843
[9,] 0.50 0.06666667 0.8571429 0.6666667 0.7500000
[10,] 0.55 0.07000000 0.8529412 0.6444444 0.7341772
[11,] 0.60 0.07666667 0.8666667 0.5777778 0.6933333
[12,] 0.65 0.09000000 0.8461538 0.4888889 0.6197183
[13,] 0.70 0.08666667 0.9130435 0.4666667 0.6176471
[14,] 0.75 0.09333333 0.9047619 0.4222222 0.5757576
[15,] 0.80 0.09333333 1.0000000 0.3777778 0.5483871
[16,] 0.85 0.10000000 1.0000000 0.3333333 0.5000000
[17,] 0.90 0.12000000 1.0000000 0.2000000 0.3333333Análisis: Para elegir el tau óptimo, se busca un equilibrio. La métrica F1, que es la media armónica de la precisión y el recall, es una excelente opción para datos desbalanceados. El tau que maximiza F1 sería una buena elección.